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ABSTRACT 



We present a detailed study of the cosmological evolution in general vector-tensor theories 
of gravity without potential terms. We consider the evolution of the vector field throughout 
the expansion history of the universe and carry out a classification of models according to the 
behavior of the vector field in each cosmological epoch. We also analyze the case in which 
the universe is dominated by the vector field, performing a complete analysis of the system 
phase map and identifying those attracting solutions which give rise to accelerated expansion. 
Moreover, we consider the evolution in a universe filled with a pressureless fluid in addition 
to the vector field and study the existence of attractors in which we can have a transition 
from matter-domination to vector-domination with accelerated expansion so that the vector 
field may play the role of dark energy. We find that the existence of solutions with late-time 
accelerated expansion is a generic prediction of vector-tensor theories and that such solutions 
typically lead to the presence of future singularities. Finally, limits from local gravity tests are 
used to get constraints on the value of the vector field at small (Solar System) scales. 
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1 Introduction 



General Relativity (GR) is considered nowadays the standard theory of the gravitational phe- 
nomena. In this theory, the gravitational interaction is mediated by a pure spin two field which 
allows to describe the geometry of the background space-time where the physical processes take 
place. In spite of its elegance and experimental success, the possible existence of other met- 
ric theories of gravity compatible with Solar System experiments as well as the importance of 
knowing the true theory of the gravitational interaction to describe properly relativistic astro- 
physical objects motivated long time ago the search for alternative gravitational theories with 
new degrees of freedom, in addition to the metric tensor. The first of these attempts was the 
so-called Brans-Dicke theory in which, apart from having the metric tensor, the gravitational 
interaction is described by a scalar field which makes the gravitational constant to vary in space 
and time. Since the time when this new theory was proposed in the early 60's, many other 
possibilities containing new degrees of freedom (not necessarily scalars) have been proposed. 

In cosmology, the interest in alternatives to GR has grown recently because of the exis- 
tence of several problems which cannot be satisfactorily addressed within the framework of GR 
without resorting to unknown components, namely, dark matter and dark energy. The need of 
new dark matter particles first appeared in order to account for the rotation curves of galaxies 
which seemed to indicate the existence of some kind of non-luminous matter in the haloes of 
such galaxies. However, its presence is also crucial to explain some cosmological observations 
as, for instance, the formation of the structures that we observe today (galaxies, clusters, etc.). 

The other dark component, dark energy, would be responsible for the accelerated expansion 
of the universe first inferred from observations of distant type la supernovae pQ and subsequently 
confirmed by more precise supernovae measurements [2] as well as other cosmological probes, 
mainly CMB temperature power spectrum and Baryon Acoustic Oscillations [3]. Although it 
is possible to explain the accelerated expansion by introducing a cosmological constant term in 
Einstein's equations, this poses a problem from the theoretical point of view because its value, 
according to observations, is extremely small compared to the natural scale of gravity set by 
Newton's constant. Although this fact would not need to be a problem, the actual situation is 
that a theory with two scales differing in many orders of magnitude does not seem very natural, 
and that is why this problem is usually referred to as the naturalness problem. For this reason, 
many models trying to play the role of dynamical dark energy have been developed during the 
last decade and can be broadly classified in two classes: on one hand, models in which dark 
energy is a new field to add to the standard composition of the universe jl] and, on the other 
hand, models in which the accelerated expansion would be an effect of a modified theory of 
gravity [5] . Nevertheless, this distinction is not always very clear because some modified gravity 
models require the introduction of new fields which, indeed, may play the role of dark energy. 

One class of alternative theories of gravity with an extra field, in addition to the metric 
tensor, are the so-called vector-tensor theories, in which the gravitational action is modified 
by adding a vector field that is non-minimally coupled to gravity. The study of vector-tensor 
theories to describe the gravitational interaction as alternatives to GR started long time ago 
with the works by Will, Nordtvedt and Hellings [B] in the early 70 's as candidates to produce 
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preferred frame effects. After these pioneering works, vector-tensor theories were abandoned 
because gravitational experiments seemed to rule out such preferred frame effects. Moreover, 
fluctuations of the vector field could be either timelike or spacelike so that those models were 
generally thought to present instabilities. A detailed treatment on this issue was developed in 
[7] for the case without potential term and in [5] for the case with a potential term and it was 
shown that there is still some room in the parameter space for stable models. A special class 
of these theories reemerged due to the increasing interest in models with Lorentz violation [9]. 
The breaking of Lorentz invariance was achieved by the presence of a vector field whose norm 
was forced to be constant by means of a Lagrange multiplier. Moreover, some of these models 
are free of instabilities as it was shown in [10J. 

More recently, after the discovery of the accelerated expansion of the universe, vector-tensor 
theories supplemented with potential terms have got much attention as possible candidates 
for dark energy [UJ. Nevertheless, vector-tensor theories even in the absence of potential 
terms have also turned out to be compelling candidates for dark energy [T2] and they could 
even solve the aforementioned naturalness problem that most dark energy models exhibit. A 
remarkable example is the own electromagnetic field. Indeed, it has been recently shown that 
electromagnetic quantum fluctuations generated during an inflationary epoch at the electroweak 
scale can give rise to an effective cosmological constant whose value is in agreement with the 
observed value [13]. Finally, vector fields have also been used for other cosmological purposes 
as candidates to drive an inflationary epoch [TJ], to generate non-singular cosmologies [15], as 
dark matter candidates [16] and as candidates to solve some of the observed anomalies in the 
CMB power spectrum [17] . 

Given the increasing interest in vector-tensor theories of gravity in order to understand the 
wide variety of cosmological problems explained above, it seems useful to provide a detailed 
analysis of the cosmological evolution in a general vector-tensor theory so that one can know 
whether a particular theory will be able to play a determined role in the universe history. In 
this work, we shall focus on the problem of dark energy and obtain the necessary conditions 
to produce late-time acceleration, although the results obtained will be general enough to be 
useful in other cosmological contexts as well. 

The paper is organized as follows: in section 2 we describe the general vector-tensor theory 
and give the corresponding equations. Section 3 is devoted to study the evolution of the vector 
field in an isotropic Friedman- Robert son- Walker universe. We shall give the behavior of the 
vector field in the different phases that the universe has undergone, namely: inflation, radiation 
domination and matter domination. After that, in section 4 we shall study the case in which 
the temporal component of the vector field dominates the universe and, in section 5, we shall 
obtain the region in the parameter space in which we might have attracting solutions with 
accelerated expansion. The case of a universe filled with matter plus the vector field will be 
studied in section 6 where it is shown that solutions exist which describe a transition from 
matter-dominance to vector-dominance with accelerated expansion. In section 7 we use current 
limits on local gravity tests to obtain constraints on the vector field at small scales. 
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2 Vector-tensor theories of gravity 



We shall start by writing the most general action for a vector-tensor theory without any other 
restriction apart from having second order linear equations of motion [IB] : 



S[g^, A^] = J d 4 xy 



16nG 



' R + uRA^A" + aR llv N i A v + tV„4, VM" + eF^F^ 



, (1) 



with to, a, t, e dimensionless parameters and F^ = d^A u — d u Af l . In the so-called iEther- 
Einstein models, the vector field norm is fixed by introducing a lagrange multiplier in the 
action of the form \ m {A^A^ ±m 2 ) so that A^ is constrained to be either timelike or spacelike. 
In other cases, the vector field is supplied with a mass term or even more complicated potential 
terms. As we said above, throughout this work we shall focus on vector-tensor theories without 
lagrange multipliers nor potential terms. Notice also that terms given in ([1]) are the only 
possibilities without introducing new scales in the theory and that gives rise to linear equations 
of motion for the vector field. If such an action is just a low energy limit of some underlying 
theory, one would expect to have corrections involving terms of dimension higher than 4 which 
would be suppressed by some scale M. 

For subsequent calculations we shall work with an alternative form of action ([1]) obtained 
via an integration by parts: 



Slg^Af] = J d 4 x^g 



16ttG 

where the new parameters relate to the old ones as follows: 



a = a — r 

A = T 

2e + r 



R + uRA^A 1 + aR^A»A v + \{V^) 2 + eF^F^ 



(2) 



(3) 



We prefer the new form of the action because it allows a more suggestive interpretation for 
each term, namely: the e-term is nothing but the U(l) gauge invariant kinetic term for the 
vector field, the A-term is analogous to the gauge fixing term introduced in the electromagnetic 
quantization and, finally, both the u and a terms are non-minimal couplings to gravity and 
play the role of effective mass terms for the vector field driven by gravity. 

The gravitational equations obtained from action (J2J) by varying with respect to the metric 
tensor can be written in the following way: 

G,* = 8nG {uT; y + aT; v + XT^ + eT^ + T» G ) , (4) 

where T^ U G is the energy-momentum tensor corresponding to other fields rather than A^ (gen- 
erally the inflaton, matter and radiation) and we have defined: 

T£ u = 2 [DA V + A2 G»u + RA„A V - V M V^ 2 ] 
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T; v = [V a Vf3 (A a AP) - R a/3 A a A /3 ] + □ (A,A U ) - 2V a V (M {A v) A a ) + AA a A {fl R v)ct 



= 9^ [{V a A a f + 2A a V a (V^)] - 4A (M V.) {V a A a ) 



rpNG 
flU 



2 SS NG 



flU 



(5) 



with □ = V^V 4 , A 2 = A^A^ and brackets in a pair of indices denoting symmetrization with 
respect to the corresponding indices. 

Apart from the gravitational equations we can obtain a set of field equations for A^ by 
varying the action with respect to the vector field to give: 



2eV v F ta/ - \V»(y u A u ) + uRA» + aR£A v = 0. 



(6) 



As we are interested in the cosmological evolution of the vector field (especially as candidate 
for dark energy) we shall focus on the simplest case in which the field is homogeneous, i.e., we 
shall study the evolution of the zero Fourier mode of the vector field. Actually, this will 
correspond to all Fourier modes whose physical wavelengths are much larger than the Hubble 
radius (super-Hubble modes). In fact, this is the relevant part of the field for the cosmological 
expansion evolution, although the inhomogeneous part could also be very important for the 
anisotropics of the CMB or for structure formation, but this is out of the scope of this work. 
Besides, we shall choose the spatial component of the field lying along the z-axis in such a 
way that we can write A^ = (Ao(t), 0, 0, A z (t)) and, therefore, we have axial symmetry around 
that axis. Thus, the metric tensor will be appropriately described by that of the axisymmetric 
Bianchi type I space-time: 



ds 2 = dt 2 - a ± (t) 2 (dx 2 + dy 2 ) - a\\(t) 2 dz 2 . 



(7) 



where a± and a\\ are the transverse and longitudinal scale factors respectively. For this metric, 
the field equations read: 



A 



A + (2H ± + H ll )A 



+ 



(2cj + a)(2Hl + HfA + (2cu + a + X)(2H ± + H\\) + 2u{2H ± H\\ + H\) A = 0, 



2e 



A Z + {2H^-H,)A Z 



+ 



(8) 



(9) 



2u(3Hl + 2H ± ) + (2u + a)(H\\ + H$ + 2H ± H\\ 



A z = 0, 



where a dot stands for derivative with respect to the cosmic time t and H\\ = ay /ay and 
H± = d_ L /a± are the longitudinal and transverse expansion rates respectively. In these equations 
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we see that the expansion of the universe provides and effective mass for each component of 
the vector field as well as a friction term. It is interesting to note that the dynamics of A$ and 
A z are driven by A and e respectively, whereas the rest of parameters of the action, u> and a, 
only affect the effective mass of the field. In fact, the presence of no n- vanishing values for A 
and e ensures the existence of evolving A and A z respectively. On the contrary, if one of these 
parameters is zero, the corresponding component does not have dynamics and, in general, will 
vanish. 

The highly isotropic CMB power spectrum that we observe today shows that the anisotropy 
at the last scattering surface was very small so that it is justified to consider small deviations 
from a pure isotropic universe so that 

a±(t) = a(t) 

a„ (f) = a(t)(l + ±h\ (10) 

with /}<1 and a(t) the isotropic scale factor. Then, we can describe the anisotropy by means 
of the degree of anisotropy h in terms of which the metric can be written as a perturbed FLRW 
given by: 

ds 2 = dt 2 - a{t) 2 (S i:j + h i:j ) dx l dx j (11) 

with hij = h5i Z 5j Z and a(t) the usual scale factor. Moreover, the degree of anisotropy can be 
related to the linearized Einstein tensor for the Bianchi type I metric as follows: 

g -- g >=M( a ) +o ^ < i2 > 

where G± = G x x = G v y and G\\ = G z z . Then, from Einstein equations we can obtain the evolution 
for the degree of anisotropy which happens to be: 



h = 16nG J J a 3 (p|| - p ± )dt 



dt (13) 



In principle, the problem has not been solved yet because the expression inside the integral 
will depend on h as well. However, we can obtain an approximate solution by replacing such 
an integrand by its expression in the isotropic case. In fact, one could obtain more accurate 
solutions by an iterative process. Thus, if we now assume that the only source of anisotropy 
comes from the spatial component of the vector field we have that: 



a 2 (P\\-P±) = (24u; e + 12u t a e + 8a t + 3a 2 )H 2 + (12u e + 6u e a t + Aa t + o 2 )H 



A 



2 



+4a t HA z A z - 2(2 + a t )A 2 z (14) 

where H = a/a the Hubble expansion rate and we have introduced the notation u t = u/e 
and a e = er/e. In the rest of this work we shall perform a detailed analysis of the isotropic 
evolution so that we could eventually use those results to evaluate ([TBI and, therefore, to 
discriminate those models in which the degree of anisotropy grows or decays as the universe ex- 
pands, that is, what models would give rise to large scale anisotropies. For a detailed treatment 
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on the anisotropy evolution in dark energy models see [19]. Moreover, the generated large-scale 
anisotropy would affect the photons coming from the last scattering surface so that it would 
give a new contribution to the low multipoles of the CMB. In fact, this can be used to rule out 
those models in which the new contribution is larger than the observed one. 

To end this section we would like to comment on a very interesting feature of these models. 
In [20] it was shown that a dark energy field carrying a non-vanishing density of momentum 
could modify the usual interpretation of the CMB dipole as well as the value of the quadrupole. 
The density of momentum of dark energy can be interpreted as a relative motion of this com- 
ponent with respect to the others. Given that photons, baryons and dark matter particles were 
strongly coupled in the early universe, they originally shared a common large scale rest frame so 
that, after becoming decoupled, they should remain at rest with respect to each other because 
of momentum conservation. However, the presence of a dark energy fluid with a non-vanishing 
density of momentum allow the existence of relative velocities among all the components with- 
out violating the momentum conservation. Therefore, as a vector field has spatial components 
one would expect it to carry density of momentum, essentially determined by A{, and, as con- 
sequence, it would be a natural candidate for a moving dark energy model. Nevertheless, once 
one uses the equations of motion it turns out that the density of momentum vanishes identi- 
cally, i.e., T° l = over the field equations. Hence, although the vector field can provide large 
scale anisotropy supported by its spatial component, it does not modify the cosmic rest frame 
so that the rest of fluids will generally share a common large scale rest frame and no effects 
on the CMB dipole are expected. However, it is possible that vector perturbations could be 
supported by perturbations of the vector field so that, unlike in standard ACDM, we could 
obtain large peculiar velocities at large scales as those observed in [2"T] . 



3 Evolution in an isotropic universe 



In this section we shall write down and solve the equations for the vector field in a universe 
dominated by an isotropic perfect fluid, with the energy density of the vector field negligible. 
In such a case, both expansion factors become the same a± — a\\ — a as well as the expansion 
rates H± = H\\ = H. This is equivalent to neglect the effects of the small degree of anisotropy 
that may be present. In such a case, the field equations read: 



A + 3HA + 



3{4u x + a x )H 2 + 3(1 + 2u x + a x )H 



A z + HA Z + - 



A 




0. 



3(Au e + a e )H 2 + (<ouo e + a e )H 
On the other hand, the energy density associated to A and A z in the isotropic case are: 

PA = A 



(15) 



3(3 + 2uo x + 2a x )H 2 A 2 + 6(1 + 2uo x + a x )HA A + A 2 
2e 



Pa z 



-(3u e + cr e )H 2 A 2 z + (6u e + a € )HA A + A 2 C 



(16) 
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In these expressions we have introduced again the notation oj\ = ou/X, a\ = a/A, u t = u/e and 
a e = a/e, which only makes sense for A ^ and e ^ 0. However, if that is not the case, the 
corresponding component does not have dynamics and, generally, vanishes so that it does not 
play any role, as we explained above. 

Although we are restricting ourselves to the case of isotropic expansion, one should be 
aware of the fact that the presence of a small shear could modify the evolution of the spatial 
components of the vector field because, in the stability analysis around a FRW metric, there 
could be vanishing eigenvalues [22J and, as a consequence, the evolution of A z could be different 
from that determined by ffl5|) . 

In next subsections we shall study the evolution of the vector field in the different epochs 
of the expansion history of the universe and carry out a classification of the models according 
to their behaviors. 



3.1 Inflationary (de Sitter) epoch 

During the inflationary era, the universe undergoes an exponential expansion so that the Hubble 
parameter H is constant, i.e., a oc e Ht . In such a case, the solutions to (IT51) can be expressed 
as: 

A = (A+e coHt + A^e- coHt )e- 3Ht / 2 

A z = {A + z e c * m + A^e- c * m )e- m l 2 (17) 



where 



c = - 48cu A - 12a A (18) 

1 



-VI - 24cj e - 6a e (19) 

Then, the evolution of the temporal component depends on whether cq is real or complex. That 
way, we obtain that it oscillates with frequency |c | and it is modulated by a damping factor of 
the form e~ 3Ht ^ 2 for 16cj a + 4a A > 3 whereas it evolves as e ( c o-3/2)m f Qr iq Ux _|_ /\ l(Jx < 3 When 
16o;a + 4cr A = 3 we have that c = and the vector field evolve as A = (C° + (7° t)e~ 3Ht ^ 2 . 

Concerning A z , it has an oscillating evolution with frequency \c z \ suppressed by e~ Ht ^ 2 for 
24u; A + 6a A > 1 and it evolves as e (Cz ~ 1/2)H * for 2Au x + 6a A < 1. Finally, for 2Au x + 6er A = 1 we 
have that c z = and the field evolve as A z = (Q + Cf t)e- 3Ht / 2 . 

When we insert solutions (fT7j) into ffTB]) we obtain the energy density evolutions, which can 
be written as: 



PAo = p At a 2 ^- 3 + p A -a- 2 ^ 3 (20) 
Pa, = p At a 2c >- 3 + p A7 a- 2c *~ 3 (21) 
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where we have defined: 



Pa ± 



A [3(3 + 2uo x + 2a x ) + 6(1 + 2uo x + a x )(±c - 3/2) + (±c - 3/2) 2 ] (HA±f (22) 
-2e [-(3c* + a e ) + (6c* + a e )(±c z - 1/2) + (±c 2 - 1/2) 2 ] (HAf) 2 . (23) 



In these expressions we see that the temporal component is suppressed during inflation in 
models with 2c < 3 =>- 4u x + a x > whereas an inflationary epoch amplifies A for those 
models with 4u x + a x < 0. For 4u x + a x = 0, the temporal component has constant energy 
density. Moreover, for the aforementioned special case with c = the energy density is given 
by: 



PA \c =0 

so it decays as ~ ta~ 3 (unless 5 — 8uj x = 0). 



^C 2 ° [3(5 - 8tu x )HC° t - (2Acu x - 15)HC° + 2C 2 °] 



(24) 



Analogously we obtain that the spatial component is amplified during inflation in models 
with 12u e + 3cr e < —4, it is suppressed for those models with 12lu € + 3a e > —4 and it has 
constant energy density if the condition 12u; e + 3cr e = —4 is satisfied. Again, for the special 
case with c z = we have a different evolution given by: 



Pa z 



3C2 [(5 



\2u t )HC z 2 t - (12u e - b)HCl + 2C 2 2 



,-3 



(25) 



so it decays as ~ ta (unless 5 — 12uj e = 0). 



Finally, notice that the oscillating behavior of the field will translate into an oscillating 
evolution of the energy density so that in those cases in which the field oscillates the energy 
density for the corresponding component is suppressed by a factor a -3 . In fact, if the inflationary 
era lasts iV e-folds, i.e., the scale factor increases as a en d/cii n = e N , we can calculate the 
amplification or suppression of the field at the end of inflation and is given by: 



In 



In 



A {t 



end I 



Ao(ti n ) 
A-z (tend) 



A z {ti 



Re(c ) - - 
Re(c,) - \ 



N 



N 



For the energy densities of each component we can proceed similarly to obtain: 



In 



In 



PAp (tend) 
PA (tin) 
PA Z (tend) 
. PA Z (tin) 



= [2Re(c ) - 3)] N 
= [2Re(c z )-3)]N. 



(26) 
(27) 

(28) 
(29) 



In these expressions, the real parts of c and c z have to be taken because they are either real 
or purely imaginary and, in the latter case, the only effect of co jZ is to produce oscillations, 
but neither suppression nor amplification of the vector field, as we commented above. In the 
special cases with c z = we have to add a term In t -f A on the right hand sides of the above 
expressions, although such a term is usually much smaller than N and can be safely neglected. 
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3.2 Barotropic fluid domination 



In a universe dominated by a barotropic perfect fluid with constant equation of state w = p/p, 
the scale factor evolves according to a power law of the form a oc t p with p = 37^3 so that 
H = p/t. In such a case, the field equations ffl~5]) have the following solutions: 



A+t a+ 



A (t) 

A z (t) = A+tP+ + A~tP 



with 



a± 



1 r 

2 
1 r 



1 - 3p ± v 7 ! + 6(1 + 4u; A + 2a x )p + 3(3 - I6cu x - Aa x )p 2 



l-p±y/l + 2(6lu € + a e - l)p - (2Au € + 6a e - l)p 2 



(30) 

(31) 
(32) 



As in the inflationary epoch, we see that the components of the vector field have two different 
types of evolution depending on whether a± and (3± are real or complex. Thus, if the term 
inside the root is positive, the corresponding component of the field will evolve as a power law 
essentially given by the growing mode whereas if the term inside the root is negative, the field 
will oscillate with an amplitude proportional to t( 1_3p )/ 2 for Aq or £( 1-p )/ 2 for A z . In fact, for 
a,/3 G C, the solutions of the vector field can be expressed as: 



A (t) = t Re ^ [C? cos (Im(a) hit) + C 2 ° cos (Im(a) hit 
A z (t) = t Re(/3) [CI cos (Im(/?) In t) + C z 2 cos (Im(/3) In t)} 



(33) 



where we see that the vector field actually oscillates harmonically in In t and not in the proper 
time t. There is still another special case which is that when we have degeneration in the solu- 
tions, i.e., when a + = a_ or (3 + = /?_. If any of these relations takes place, the corresponding 
component of the vector field has a logarithmic solution in addition to the potential solution, 
being the complete solution as follows: 



A 

A, 



(C 1 ° + C 2 °lnt) t^/ 2 
{CI + C 2 2 lnt) t {1 - p)l2 



(34) 



The evolution of the energy densities are achieved by inserting solutions (130!) into (|T6|) and, 
for a, (3 6 M, are given by: 



where: 



PAt 



PA =P A +a K « +p A -a K o 
Pa z = P A +a Kt + p A -a K7 



A [3(3 + 2u x + 2a x )p 2 + 6(1 + 2uo x + a x )p a ± + a 2 ± ] (A±) 2 
-2e [-(3w e + a e )p 2 + (Qcu e + a e )p (3 ± + (3 2 ± ] (Af) 2 



(35) 
(36) 

(37) 
(38) 
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and 



4 = 2^ 



= 2 &Vp 

P 

Again, given that Kq z > k,q z , the energy densities will evolve proportionally to a K °> z so that the 
vector field behaves as the superposition of two perfect fluids with equations of state u> ,2 = 
— |(1+Kq^). However, in some cases, it might happen that p A + vanishes and, as a consequence, 

the corresponding energy density will evolve as a K °- z rather than a K ^ z . 

When a± and j3± are complex, the energy densities still evolve as p& >z = Pa+ aK °' z although 
we must replace a + — > Re(a+) and (3+ — > Re(/3+) in (1391) and p A + are oscillating functions 
instead of constants. 

For the degenerate case with a logarithmic solution, the energy densities are given by: 
(C° 2 ) 2 



PA()\a-\-=a— ^ 
PA z \p+=l3^ = 



2(1 -p) 

'2 



(cp 2 



l-3p) 



[l-2p+(24u x -l5)p 2 ] f^ + \nt\ + 2(1 - p) 
[-1 +6p+ (12u e - 5)p 2 ] (j± +\ntj+ 2(3p - 1) 



a K0 (40) 
a K * (41) 



Therefore, the evolution is the same as in the non-degenerate case modified by a logarithmic 
variation. 

In any case, we can carry out a classification of the models according to whether the energy 
density of each component grows, decays or remains constant. Finally, we can also identify 
scaling behaviors or whether the energy density of the vector field grows or decays with respect 
to that of the dominant component. In Fig. [Hand Fig. [2] we show the evolution in the different 
regions in the parameter space. 



3.2.1 Radiation dominated epoch 

In the radiation dominated epoch we have that p = 1/2 and the evolution of the vector field, 
according to fl32|) . is given by: 



a 



R 1 

± 4 



— [1 =F V25 + \2a x \ (42) 



/?£ = i[i±vT=2^ (43) 

These expressions allow to obtain that the temporal component evolves as a power law for <j\ > 
— 1| whereas it oscillates with an amplitude decaying as t -1 / 4 for a\ < — 1|. The degenerate 
case happens for a x = — ft. For the spatial component of the vector field, we have power law 
behavior for a t < \ and it oscillates with an amplitude growing as t 1 / 4 for a e > \, whereas 
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the logarithmic solution appears for a e = |. Notice that the evolution of the vector field does 
not depend on the parameter u of the action, which is due to the fact that the Ricci scalar R 
vanishes in a universe dominated by radiation. 

Concerning the evolution of the energy densities, according to (|39|) . we have that: 

re£t = -5 ± V25 + 12a x (44) 
k* ± = -5 ± VI - 2a e (45) 

If we want the vector field to be dominated by its temporal component we must impose the 
condition k 0+ > k z+ , which leads to the constraint 6a\ + a e > —12. That way we prevent the 
generation of large scale anisotropy that would be in conflict with observations. 

The temporal component of the vector field will have constant energy density in this epoch 
for k r + = which is satisfied by the model with a\ = and it has scaling evolution if k r + = —4 
which happens for a\ = —2. From this last condition we also obtain that, in models with 
<j\ > —2, the energy density associated to the temporal component grows with respect to that 
of the radiation fluid pn whereas in those models with a\ < —2, the ratio Pa /pr decays as the 
universe expands. See Fig. [1] for a summary of these behaviors. 

Concerning the spatial component, the condition of constant energy density is satisfied for 
models with a t = —12 whereas it scales as radiation in models with a e = 0. Finally, the energy 
density associated to the spatial component of the vector field grows (decays) with respect to 
Pr in models with a t < (cr e > 0). This classification is shown in Fig. [2j 



3.2.2 Matter dominated epoch 

In a universe dominated by a pressureless fluid the scale factor evolves as a oc t 2/3 so that 

= ~ [3 T V81 - A8u x + 24a A J (46) 

/3f = \ [I ± VI - 24cu e - 12a e j (47) 
o 

In this case, the temporal component oscillates with amplitude proportional to i -1 / 2 for 1Quj\ — 
8a\ > 27, it evolves with a power law for 16lj\ — 8a\ < 27 and the degenerate case corresponds 
to 16u)\ — 8a x = 27. On the other hand, A z follows a power law evolution for 24uj e + 12cr e < 1, 
it oscillates with amplitude proportional to t 1 ' 6 for 2Au e + 12a e > 1 and the degenerate case 
happens for 24co> £ + 12a,; = 1. 

Moreover, k r ± and k^ ± become: 

9 ± V81 - 48cj a + 24a A J (48) 
9 ± VI - 24w e - 12a e J (49) 



k m - -\- 
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In this epoch, the condition for the temporal contribution to dominate over the spatial one 
reads Q(u e - 2u x ) + 3(a e - 2a x ) > -80. 

The scaling behavior for the energy density associated to the temporal component in this 
epoch is obtained from the condition = —3, whose solution is 2u\ — a x = 3. Moreover, for 
2u; A — cr A < 3 the energy density of the vector field grows with respect to that of matter and, 
as a consequence, the universe is eventually dominated by it. Notice that this is a necessary 
condition to have a dark energy model, i.e., those models with 2uj\ — o~ A > 3 will never dominate 



Property 


Inflation 


Radiation 


Matter 


Oscillating 


4u A + a x > f 


a A <-§ 


2cj a - a x > f 


Decaying 


< Auj x + a x < f 


-§ < vx < -2 


3 < 2uj x - <r\ < f 


Scaling 


4cu A + a x = 


a x = -2 


2uj x -a x = 3 


Growing 


4w A + a x < 


<?\ > -2 


2uj x - a x <3 


Constant 


4cu A + a x = 


a x = 


2io x -a x = 



Inflation Radiation Matter 




Figure 1: Behavior of the temporal component of the vector field in the different phases of the universe 
history according to the values of the parameters of the model. The labels Decaying and Growing 
refers to the behavior with respect to the dominant component, whereas Scaling means that it evolves 
in the same way as the background. 
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the energy content of the universe if there is a matter component and, as a consequence, it 
cannot be the responsible for the present acceleration. Fig. [1] shows this classification. 

For the spatial component we obtain constant energy density in models with 2u e + a t = 
—20/3 and scaling evolution for 2u e + a e — —2/3. Finally, the energy density of A z grows 
(decays) with respect to pu in models with 2u e + a e < —2/3 (2a> e + cr e > —2/3). This is shown 
in Fig. H 



Property 


inflation 


Radiation 


Matter 


Oscillating 


4w e + cr e > i 


a t > \ 


1 

2w e + cr e > j2 


Decaying 


-f < Ato e + a e < \ 


< a e < 1 


-| < 2iO t + (T e < ± 


Scaling 


4u; e + cr e = — | 


a e = 


2uj t + a € = - 1 


Growing 


4w e + a e < -| 


a £ < 


2w £ +cr e < -| 


Constant 


4w e + cr e = — | 


a e = -12 


2w e + a e = -f 



Inflation Radiation Matter 




Figure 2: Behavior of the spatial component of the vector field in the different phases of the universe 
history according to the values of the parameters of the model. The labels Decaying and Growing 
refers to the behavior with respect to the dominant component, whereas Scaling means that it evolves 
in the same way as the background. 
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4 Vector dominance 



In this section we shall study the case in which the universe becomes dominated by the temporal 
component of the vector field so that the anisotropy is small and we can use the isotropic 
equations. Thus, we have the field equation for Aq given in ffl5l) and the two Einstein equations: 

3H 2 = 8nGp Ao (50) 
3H 2 + 2H = -8nGp Ao . (51) 

Although, of course, only two of the three equations are independent, it will be useful to work 
with all of them. 

For the subsequent analysis, it will be convenient to introduce the field variable x = rf J^° 
so that we can obtain the following autonomous system: 



dH g (2cu A + a x )x 2 - 2(4u x + a x )x - 3(4u x + cx x )(2uj x ±^±h_ R ^ 
dN ~ ' [ x + 3{2u x + a x + l)] 2 

dx x 2 + Q(2u x + a x + l)x + 3{2u x + 2a x + 3) 



dN x + 3(2u x + a x + 1) 



x (53) 



with N = In a. These two equations can be combined to give the following equation for the 
trajectories in the phase map: 

dB_ = (2lo x + a x )x 2 - 2(Auo x + a x )x - 3(Auj x + a x ) (2uj x + a x + 1) H 

dx [x 2 + Q(2uj x + a x + l)x + 3(2uj x + 2a x + 3)}[x + 3(2LO X + a x + l)} x 1 ' 

This equation can be readily integrated for given values of the parameters, although we shall 
not do it, but we shall study the phase map and, from its features, we shall obtain the relevant 
information. 

In addition to the equations given above, we have the following constraint provided by the 
Friedman equation: 

\\A 2 [x 2 + 6(2cu A + a x + l)x + 3{2uj x + 2a x + 3)] = 1 (55) 

This relation constrains the possible values of x because the condition 

A [x 2 + 6(1 + a x + 2uj x )x + 3(3 + 2<r A + 2uj x )] > (56) 

must hold. This condition restricts the physically admissible trajectories as those for which the 
vector field carries positive energy density and can always be achieved by means of a suitable 
choice of the sign of A. Notice that all the dependency of the problem on the parameter A is 
indeed contained in this condition. 

In (|53|) we see that the equation ^ = has two solutions: H = and x = xf, with: 
± _ 4u; A + a x ± y/(4cu A + qt a ) [2{huj x + 2a x ) + 3(2u x + q A ) 2 ] 
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whereas the solutions of = are x = and x = xf, with: 

xt = -3(2u x + a x + 1) ± v/6(5^ A + 2cr x ) + 9(2u x + a x ) 2 . (58) 

Therefore, the autonomous system has generally three critical points: Pq = (0, 0) and P± = 
(xf,0). Moreover, when the equality xf = xf takes place, we have a critical line instead of 
a critical point because both ^ and vanish regardless the value of H. Notice that the 
critical points P± only exist when the constraint 6(5u x + 2cr\) + 9(2u x + a x ) 2 > is satisfied, 
which corresponds to the white region shown in Fig. [3j Moreover, the trajectories in the phase 
map have vertical tangents in the lines x = and x = xf so that, apart from being critical 
points, they are vertical separatrices (see Fig. H]). However, these are not the only separatrices 
in the phase map, but we have another vertical one in x s = —3(2u x + a x + 1) as well as one 
horizontal one in the horizontal axis H = 0. One interesting feature of the phase map is that 
x s = (x+ + x~)/2, i.e., the separatrix at x s is always located in the middle of the two critical 
points. Notice that the critical points xf also separate the region with physically admissible 
trajectories according to fl56|) which imposes the energy density of the vector field to be positive, 
although to identify each of these regions we need to specify the sign of A. According to the 
previous discussion, the phase map will be divided in several rectangular regions parallel to the 
axis which, indeed, are disconnected from each other, i.e., the trajectories will not be able to 
cross from one to another. However, the particular picture will depend on the particular values 
of the parameters, being possible to distinguish the following cases (see Fig. Hj): 

• Case I: The critical points xf exist and are different from each other, which imposes the 
condition 9(2u x + a x + l) 2 — ?>{2uj x + 2a x + 3) > and corresponds to the grey region in 
Fig. [3j Moreover, in this case we still have three different possibilities: 

— Case la: Both critical points are different from zero and they do not coincide with 
x s . In this case we have 4 vertical separatrices and the phase map is divided into 10 
disconnected regions. 

— Case lb: One of the two critical points is zero. This case happens when 2uj x + 2a x + 
3 = but 2uj x + a x + 1 ^ so that they are not zero simultaneously. In fact, if 
2u x + a x + 1 is positive (negative), then x^ (x~) is at the origin. Notice that this 
simply says that, given that the separatrix located at x s is in the middle of X- and 
x+, the critical point that is at the origin depends on the sign of x s . In this case 
we only have 3 vertical asymptotes and, as a consequence, the phase map contains 
8 regions. This case corresponds to the dashed line in Fig. [31 which represents the 
equation 2u x + 2a x + 3 = 

— Case Ic: The vertical separatrix satisfies x s = so that the two critical points are 
symmetric with respect to the origin. The condition x s = reads in terms of the 
parameters 2u x + a x + 1 = and, as a consequence, the critical point Pq is absent. 
In this case we also have 3 vertical asymptotes and 8 regions in the phase map and 
corresponds to the dotted line inside the grey region in Fig. |3j 
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• Case II: Neither of the two critical points xf exists. This case satisfies 9(2u>\ + a x + l) 2 — 
3(2uj\ + 2a \ + 3) > and corresponds to the white region in Fig. [31 Again, we have 
several subcases: 

— Case Ha: The separatrix x s is different from zero. Then, we have two different 
vertical asymptotes and the phase map contains 6 regions. 

— Case lib: The separatrix is located at the origin, x s = so that 2u x + a x + 1 = (Po 
is not a critical point) and we are left with just one vertical asymptote that divides 
the phase map in 4 regions. This case corresponds to the dotted line inside the white 
region in Fig. [3j 

• Case III: The two critical points become equal so we only have one critical point which, 
indeed, coincides with 3(2u\ + a x + 1). This case is represented 
by the solid line in Fig. [3J The two possibilities we have in this case are: 

— Case Ilia: The critical point is different from zero so that we have two vertical 
asymptotes and the phase map contains 6 regions. 

— Case Illb: The critical point is zero. Then, we only have one vertical separatrix and 
four regions. This case corresponds to the particular model (a x = —2,uj x = 1/2) 
and is represented by the orange dot in Fig. [3j 

Now that we know the arrangement of the phase map for the different models attending to 
the existence and location of the critical points, we shall study the particular features of each 
critical point according to the values of the parameters: 

• P = (0,0). The eigenvalues for this critical point are: 



Thus, we have that the critical point is a saddle point for models whose parameters are 
between the lines 4u x + a x = and 2u x + 2<T\ + 3 = whereas it is an attractor node if 
the parameters are in the external region. Notice that this critical point does not exist 
when 2oj\ + a x + 1 = 0, that corresponds to Case Ic in which the separatrix x s is placed 
at the origin. When Auj x + a x = we have that fin — and the critical point becomes a 
critical line because x = is a singular point irrespectively of the value of H. 

• P± = (xf, 0). For these critical points the eigenvalues can be expressed as /i# = — (2xf + 
3) and \i x = —2xf. Then, if the critical point is positive we have an attractor node whereas 
it behaves as a repelling node if xf < —3/2. Finally, in the range —3/2 < xf < we get 
a saddle point. These ranges correspond to the regions in the parameter space showed in 
Fig. [5j The eigenvalues vanish when the corresponding critical point is located at the 
origin, i.e., 2<j\ + 2u\ + 3 = 0. Then, we can have three different cases depending on the 
sign of 2u\ + (T\ + 1, namely: 



4o; A + ct a 



2cu x + 2o- A + 3 
2uj x + a x + 1 



(59) 



2u x + a x + 1 ' 
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i) 2oj\ + o\ + 1 > 0. In this case //+ = and [i x = — 6(2u\ + <j\ + 1) < 0. 

ii) 2uj\ + a\ + l = 0. This is the case Illb described above in which both critical points 
are located at the origin and, as a consequence, /x+ = \i~ = 0. 

in) 2uj\ + &\ + 1 > 0. This case is the opposite to i), i.e., \i~ = and /i+ = — 6(2u\ + 
a x + 1) > 0. 




-4 -2 2 4 



Figure 3: This plot shows the different cases explained in the main text in the parameter space. 
The green shaded region corresponds to models in which both x± exist (Case I) whereas in the white 
region neither of them is present (Case II). The solid line identifies the models for which x+ = x_ 
corresponding to Case III. The dashed line (whose equation is 2u\ + 2o\ + 3 = 0) represents those 
models that have either x_ or x+ at the origin (Case lb) and the models whose parameters lie on 
the dotted line (with equation 1lo\ + o\ + 1 = 0) have the separatrix x s at the origin. Thus, in the 
region above (below) that line x s is negative (positive). Case Ic corresponds to the piece of the dotted 
line inside the white region whereas the piece of the dotted line inside the grey region corresponds to 
Case lib. Finally, the orange dot with parameters u>\ = 1/2, ax = —2 gives the Case Illb in which 
x± = x s = 0. 
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On the other hand, the eigenvalues fi H vanish in the case that the critical points are such 
that xf = —3/2 which is satisfied for models with a x + 4u x = 3/4. In this case we also 
obtain three possibilities in terms of the sign of 12uo x — 15/2 as follows: 

i) 12uo x — 15/2 > 0. In this case we have /% = and = 2Auo x — 15 > 0. 
it) 12cg>a — 15/2 = 0. In this case we have = /i^ = 0. 
in) 12uj\ — 15/2 < 0. In this case we have = and = 2Au x — 15 < 0. 

Separatrix x = x s . For values of x close to x s , i.e., x = x s + 5 with 5^0 the equations 
become: 

dH H 

m " "'p (60) 

with 

fi a = 9(2lu x + a x + 1) [12^ + a x (A + 3a x ) + 2lo x (5 + 6a x )] . (62) 

Then, the separatrix will attract the trajectories of the phase map for models in which 
fi s > whereas for models with fi s < the trajectories will go away from x s . Note that 
this is true for both sides of the separatrix. 



H o,o- :-— 





Figure 4: These plots shows two examples of phase maps for the autonomous system describing a 
universe dominated by the vector field. The left panel corresponds to a model lying in the grey region 
in Fig. [3] in which the two critical points P± are present whereas the right panel shows the phase map 
for a model in which these two critical points do not exist. We can see that all the critical points as 
well as the separatrix are vertical tangents and the x-axis is a horizontal separatrix, as explained in 
the main text. 
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• x — > ±oo. In this case we have that ~ — x 2 so that the trajectories will always 
approach from +00 and will move away to —00, i.e., the region with large values of x 
+00 repel the trajectories of the phase map whereas they become attracted by the region 
with x — > —00 irrespectively of the values of the parameters in the action. This means 
that the region with large positive values of x is always unstable for any choice of u\ and 
<j\ and, on the contrary, the region with negative large values of x is always stable. 




Figure 5: In this plot we show the different regions obtained in the parameter space according to 
the features of the critical points. We have shaded (in orange) the regions where the eigenvalue of 
x is negative so that the trajectories in the phase map approach the corresponding point. We have 
also indicated the sign of each eigenvalue in the form (/j, x ,ij,h). For the separatrix we indicate the 
regions where the trajectories approach the separatrix. Finally, the two plots on the right show the 
13 regions explained in the text (and summarized in table 1). The last plot corresponds to a zoom of 
the blue-shaded region in the third panel. 
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Region 


Po 




Separatrix 


P+ 


I 


Attractor node 


Saddle point ^ > 


Repelling 


Attractor node 


II 


Saddle point \i x > 


Attractor node 


Repelling 


Attractor node 


III 


Saddle point \i x > 


- 


Attractor 


- 


IV 


Saddle point /i x < 


- 


Repelling 


- 


V 


Saddle point jjb x < 


Repelling node 


Attractor 


Repelling node 


VI 


Attractor node 


Repelling node 


Attractor 


Repelling node 


VII 


Attractor node 


Repelling node 


Attractor 


Saddle point ji x > 


VIII 


Saddle point \i x > 


Repelling node 


Attracting 


Attractor node 


IX 


Saddle point /i x < 


Repelling node 


Repelling 


Attractor node 


X 


Saddle point jj, x < 


Saddle point ji x > 


Repelling 


Attractor node 


XI 


Saddle point \i x > 


Saddle point /i x > 


Attractor 


Attractor node 


XII 


Attractor node 


Saddle point [i x > 


Attractor 


Saddle point \x x > 


XIII 


Attractor node 




Repelling 





Table 1: In this table we summarize the features of the phase map for the different regions shown in 
Fig. When a given critical point is a saddle point we give the sign of the eigenvalue corresponding 
to x so that we can know whether the trajectories approach the critical point (negative eigenvalue) or 
move away from it (positive eigenvalue) along the x-direction. 
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4.1 Accelerating solutions 



In this section we shall enumerate the necessary conditions for a vector-tensor model to lead to 
accelerating solutions. To that end we shall express the equation of state in terms of the field 
variable x: 

{4u x + 2a x + l)x 2 - 2{2u x -a x - 3)x - 3(2cu A + a x + l){2u x - a x - 3) 

W = s (06) 

[x + 3(2co x + a x + l)] 2 

Notice that this equation of state only depends on x and not on the Hubble parameter H . 
The models with accelerating solutions will be those in which w evolves towards w < —1/3. 
Moreover, as this work is intended to find models in which the vector field could play the role 
of dark energy, we shall demand that the accelerated phase is an attractor. To that end, we 
shall look at the equation of state for all the possible attracting places in the phase map as well 
as in the repelling ones. Notice that, as the equation of state does not depend on H but only 
on x we only need to require attractor or repelling properties with respect to x. For instance, 
if a critical point is a saddle point but with the trajectories going towards x — > x cr u (fi x < 0) it 
will be considered as an attractor and the opposite for a repelling point. Now, we shall study 
the existence of accelerating regimes in the phase map: 

• Pq. As we pointed out above, this critical point will be an attractor in models with /i x < 0. 
On the other hand, the equation of state for this critical point is 

2u x -a x -6 

w «< = 3 (a* + ,» + !) (64) 

From this expression we see that Pq corresponds to an accelerated phase (wp < —1/3) if 
the following condition holds: 

2 ^" 1 < (65) 
2u x + a x + l K J 

which is satisfied for models in which either 1 < 2uj x < — (1 + <j x ) or —(1 + a x ) < 2uj x < 1. 
The corresponding region is shown in Fig. [6] and, in that figure, we see that there exists a 
region in which Pq is an attractor and gives rise to accelerated expansion simultaneously. 

• P±. For these critical points, the equation of state becomes: 



w P± = - (8to x + Aa x + 3) ± - y/6(5u x + 2a x ) + 9(2o; A + a x ) 2 (66) 

Notice that if the condition of accelerated expansion is satisfied for P + , then, it is also 
satisfied for P_ since wp + > wp_. The models that lead to accelerated expansion are 
shown in Fig. [6j However, in the same figure, we see that neither of the critical points 
P± behaves as an attractor in the region where we get accelerated expansion. 

Separatrix x = x s . In this case the equation of state will evolve either to +oo or to — oo 
as the trajectory approaches the separatrix. The interesting case here is when w — > — oo 
so that we get acceleration. For x = x s + 5 with 5 — > 0, the equation of state becomes: 

w s ~ (67) 
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In Fig. [6] we show the region in which the equation of state goes to — oo, corresponding 
to the condition fi s > 0. Notice that this condition also guarantees that the separatrix 
behaves as an attractor so all the cases in which the sepratrix is an attractor give rise to 
accelerated expansion. 

• x — > ±00. As we showed above, x = —00 is always an attractor whereas x = +00 is 
always a repelling point. Moreover, as x approaches ±00 the equation of state is given 
by: 

w(x^±oo) = -(4uo x + 2a x + l) (68) 

so that it gives rise to accelerated expansion in models with 2u x + a x > —1/3. This case 
is interesting for the vector field to drive an inflationary era because it could start with 
a large value of x and, as x = 00 repels the trajectories, it would be forced to evolve 
towards smaller values of x until it reaches either P + or P . Moreover, such an evolution 
can lead to accelerated expansion as long as the condition 2u x + a x > —1/3 holds. 



After having obtained the general conditions necessary to have accelerated expansion, we 
shall study the particular solutions in which the scale factor evolves as a power of the cosmic 
time, i.e., a cx t p and we shall get some analytical solutions. For this expansion law, the vector 
field evolves according to (13"U|) . although, in this case, the parameter p must be determined from 
Einstein equations. When we introduce fl30|) in (150]) we obtain that the vector field must take 
a constant value given by: 

A°° = ^ (69) 

° v/87rGA(3 + 2u x + 2a x ) 

Notice that this value only makes sense for A (3 + 2u x + 2o x ) > 0, which can always be fulfilled 
by a suitable choice of the parameter A. Moreover, we also obtain that p is given by: 

1 + 2u x + a x {7n , 

V = — : (70) 

Auo x + a x 

With this expression we can calculate the deceleration parameter: 



aa 2uj x - 1 

goo = -— = 7^ : — r 71 

a 2 2uj x + a x + 1 



This deceleration parameter must be negative in order to have an accelerated expansion, but this 
is indeed the same condition that we found above in (|65|) when studying the critical point P$. 
This was expected because the Hubble expansion rate goes to zero for a power law expansion 
and this imposes that Aq has to be constant which means that x = 0. Thus, we have just 
obtained nothing but the analytical solutions for the trajectories approaching P along the H- 
axis. Notice that this is, indeed, the only critical point attracting trajectories with accelerated 
expansion and having a finite value for the vector field equation of state. 
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Figure 6: In this plot we show the regions in which we get accelerated solutions in a universe dominated 
by the vector field for each critical point. The darkest shaded regions are those in which the accelerated 
solutions are attractors. We can see that neither for P + nor for P_ we can obtain attracting accelerated 
solutions. 



5 Transition from matter domination to vector domina- 
tion: late-time accelerated solutions 

In previous sections we have studied a universe completely dominated by the vector field and 
obtained the necessary conditions to have accelerated solutions. In this section, however, we 
shall consider the case in which the universe contains matter in addition to the vector field and 
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study the circumstances under which we can get a transition from a matter dominated universe 
to an accelerated phase provided by the vector field so that the vector field can play the role 
of dark energy. To do so, we shall proceed as in previous sections, i.e., we shall obtain the 
corresponding autonomous system and identify attracting solutions in which the vector field 
eventually dominates the energy content of the universe and has equation of state smaller that 
-1/3. 

Before going on with the study of the autonomous system, we remind that a necessary 
condition to have a candidate to dark energy is that the energy density associated to the vector 
field decays slower than that of a pressureless fluid in the matter dominated epoch. This 
requirement guarantees the dominance of the vector field at late times so that it can drive 
the expansion of the universe. According to Table 1, such models are those whose parameters 
satisfy the condition 

2u x - a x < 3. (72) 

We shall take this condition as a necessary requirement for the model to be able to play the 
role of dark energy. 

The system of equations must be modified by introducing the matter contribution to the 
Friedman equation so that: 

3H 2 = 8nG(p Ao +p M ). (73) 
Moreover, we have a new equation provided by the energy conservation of the matter fluid: 

p M + SHp M = (74) 

As we did in the previous sections, we shall introduce the field variable x = d J" n ^° and the 
matter energy density will be described by the density parameter VL M = |^ . In terms of these 
variables, we can obtain the following autonomous system: 

dH 3 (2u A + a x )x 2 - 2(Au x + a x )x - 3(4cu A + ct\)(2u x + + 1) + F H £l M R ,^ 
dN ~ ° [ x + 3(2cu A + a x + l)] 2 - [9(2lo x + a A ) 2 + 6(5cu A + 2a x )} Vl M 
dx [x 2 + 6(2cu x + a x + l)x + 3(2cu A + 2a A + 3)] [x(x + 3(2u x + a x + 1)) + F x tt 



Ml 



dN [ x + 3(2w A + a x + 1)]' - [9(2cu x + a x ) 2 + 6(5u x + 2a x )} n 



M 



M 



3 [(4io A + 2a x + l)x 2 - (Auj x - 2a x - 6)x - 3(2io A + a x + l)(2u x -a x - 3)] 



-; — — — — — II — f2 IVf ) f2 A//- 

dN [ x + 3(2u x + cr A + l)] 2 - [9(2o; A + a x ) 2 + 6(5u x + 2a x )\ Q M 

where we have defined: 

F H = ~\ [(4u) X + 2(x x + l)x 2 -2(2u x -a x - 3)^ + 6(4^ + ^(2^ + ^+9(2^-1)] 

F x = -^[{4cu x + 2a x + l)x-2uj x + a x + 3}. (76) 

These equations are supplemented by the following constraint provided by the Friedman equa- 
tion: 

-\A 2 [x 2 + 6(2u x + a x + l)x + 3(2cj a + 2a A + 3)] + Q M = 1 (77) 
o 
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As before, this relation will determine the sign of the parameter A in order to fulfill the condition 
A [x 2 + 6(2lu x + a x + l)x + 3(2u x + 2<x A + 3)] > 0. 

The equation of state for the vector field is given in this case by: 

(Au x + 2a x + l)x 2 - 2(2lu x - a x - 3)x - 3(2lu x + a x + l)(2w A - <t x ~ 3) ,_ Q , 

w = o (78) 

[x + 3(2cu x + a x + l)] 2 - [9(2c^ A + a x ) 2 + 6(5u x + 2a x )} Q M 



The sections {x, Qm} of the phase map do not depend on H. Notice that, as expected, the 
equations for ^ and reduce to (15"3"1) when Q M — so that all the critical points analyzed in 
the previous sections are also critical points of (1751) with Qm = 0, which is the interesting case 
here because this means that the vector fields eventually dominates the energy of the universe. 
Apart from these critical points with Qm = we also have critical points with Qm = 1 which 
correspond to situations where matter drives the universe expansion. As we are interested in 
obtaining solutions leading to late-time accelerated expansion driven by the vector field, we 
shall study just the critical points with Q M = so that the critical values of x and H are 
the same as those studied in the vector dominance case. However, the features of the critical 
points may change because of the presence of matter. Therefore, for each critical point, we 
shall identify the necessary conditions for the corresponding critical point to be an attractor 
with respect to Qm and x and, then, compute the equation of state in the critical point. This is 
possible because the equations for Q M and x does not depend on H. Notice that the condition 
for Qm = to be an attractor will, in general, differ from the condition given in (1721) obtained 
by imposing that the energy density of the vector field grows with respect to that of matter. 
This is so because to achieve that condition we assumed that the amount of matter was initially 
dominant with respect to that of the vector field and, however, for Q M = to be an attractor, 
such a condition does not need to be satisfied. Finally, it is interesting to remark that the 
region with Qm > is disconnected from that with Qm < because the trajectories are always 
tangent to the plane Qm — which ensures that the energy density of matter remains positive 
as long as it is initially positive. See Fig. [7J to have an idea of how the phase maps look like. 

Let us analyze then each critical point for this case: 



• P = (0,0). For this critical point, the linearized system becomes: 

dH 

dN 
dx 
dN 
dQM 
dN 

The eigenvalues for this system are the same as those of the vector dominance case plus 
fMQ M = {2oj x — a x — 3)/{2uj x + (7 X + 1), which determines the stability of the solutions 
with Qm — ► 0. Therefore, the analysis proceeds exactly the same as before with the 
supplementary condition fiQ M < ensuring the late-time domination of the vector field. 



4u x + a x 



-H 



2u x + a x + 1 
2lu x + 2a A + 3 
2lu x + a x + l 
2u x - cr A - 3 
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2u x - cr A - 3 
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Notice that, as commented above, this condition is not the same as that given in (1721) . 
Indeed, this supplementary condition happens not to reduce the region of the parameter 
space in which we get attracting solutions with accelerated expansion, i.e., all the models 
indicated in Fig. (J5J) corresponding to these type of solutions have indeed Qm = as an 
attractor. 

• P± = (x±, 0). In this case the linearized system becomes diagonal with the same eigenval- 
ues for x and H as in the vector dominance case, i.e., /i# = — (2xf + 3) and \i x = —2xf. 
Moreover, the eigenvalue for fl M is given by /zq m = 4xf + 3. Hence, as in the case of 
Po the stability analysis is the same as that already performed above, although we must 
impose the condition 4xf + 3 < so that the vector field eventually dominates. However, 
when we impose the latter condition we find that these critical points happen not to be 
attractors for any value of the parameters. 



So far, we have seen that the presence of a matter fluid only affects the solutions in the 
sense that the parameter space is restricted to that region in which Qm = is an attractor, 
otherwise the vector field would never dominate and we cannot produce late-time acceleration. 
In other words, the features of the critical points remain the same as those studied in the vector 
domination case, although only the cases in the allowed region are admissible. The novelties 
appear when studying the separatrix: 



1,0- 



0,5 -„ 



/I 
</ i i 
"4 , 
~* ,--/ / t 



a M o,o- 



-0,5- 



-1,0- 




A \ \ 

— I 
-*^AM [ 

MMI t 

-VKM 1 t 

1 t 

-■WW 1 t 

1 1 t 



^ Wv 
VA\ 
k \ \ 

\ IWW- 
\ 1 WW*. 
\ \ 

\\\ \\\ 
V I \ \ \\v 

IS \ VV NV 

iuuw 
1 k i. \ >. \ \ 

I \ w / 
' tl II 

II M \ 



1 \ \~. 
\ \ w 
1 i \ 
1 ^ \n 
\ \ \^ 

I 1 \v 

I I \n 
1 1 \x 
1 1 w 

j j \ \ 
\{\\ 

1 UN 
1 1 \ V 

\ \ \\ 

j u\ 



r^sSs I 

J I 
1 1 

14 41 
1411 

;/ 4 11 
//J4 4 t 1 1 

4/4/11 

- — s ss s / / / / / / ; 1 1 1 

-y V J J / v 4 / / / / ft 1 t 

^/«//««/ /ftrt\ 
--*■////// j // / r r 1 1 1 

/ J 4 4 S 4 4 4 4 / / / t 1 t 

' ! 

Wlt44i//d?//t/tt1 

•—■\ \ U 

— — -N \ 4S 

^-/l- 4 I J 
-<-ss/4 I J 
Wit I J I 

-^s/4 4 I ) I 
-^s/4 4 I I I 
till) 
-'■si Hill 

~-^s4 4 4 I I I 
'/I 14 111 
-'4 4 I I It I 
— V/ t It I I 
-'■4/ I I II I 
-'Jit I II I 
4 I I I I I I 
4 I I I I I 
4 till! 
'J It I III 
i^JJ-l I 'J 



\\\\\\\\\\ 
\ \\\\\\\ 
tl 

Hi! 
iii 
j 1 1 1 1 1 



jiu 
mi 

iitt . 
jiu i 
j 1 1 1 1 1 

I j 
j 1 1 1 j i 

LUJ_U 




: n M 0,0- 



-i,o4 



z y ^. — . — 

I f /s^^ — 

t 1W. , 

t / t rr**"-~ 
t / 1 1 ft* — - 
t ( ////y^-~ 
t t f Wf-- 



I I t //w 
Will! rs"~ 
\\1\tt rrss^ 
5 ) J t It t 
\\\\ t It"- 

\\\ \ \t I r/---* 

\ U T 7 / / /y^-< 

Wit If-* 

*A\ \ I IS*^ 
.'-V f y^"^ 




i-U-P- 









'// ttt\\\ 


/ / 


1 1 1 1 ft t 


/ / 


If 


/ / 


ttttlit 




1 1 11 1 1 1 1 




11111111 




11111111 




11111111 




11111111 




11111111 
11111111 




1 1 1 1 1 1 1 1 
11111111 




11111111 




11111111 
11111111 




11111111 




11111111 








•*•> 1 

■~s.\.\^\, \ 1 
I ( 

^.sw \ I I / 
WW i I 4 
\ \ 1 1 J 4 4 
U I I 4 4 
\ \\ I I tdJ 
1 1 J //// 

i J IW 
^\\\\ \ I I 44SSS 
WW 1 1 / 4/S/f' 

s\\ 1 (//// 

— *~" 



\ ^^ 
\\\ 
\ \ \\ 

\ \ \wcv 
\ \ \\\\\ 
\ \ \\ 

\ \ \\ \ \N 
\ \ \ \ 1 \\ 
\ \ \ W \\\ 
\ \ \ \ \ \\\ 
\\\\\\\\ 
\\\\\\\\ 
\\\\\\\\ 
\ \ \ \ \ \ \ \ 
\\\\\\W 








Figure 7: In these two plots we show two examples of phase maps corresponding to the cases when 
the two critical points are present (left) and when they are not (right), or, equivalently, when the 
separatrix is open from above or from below. 
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Separatrix. In this case the separatrix is no longer given by x = x s , but by the parabola: 

O [x + 3(2^ a + <xa + 1)] 2 

UM [9(2^ + ^ + 6(5^ + 2^)] 

Notice that the vertex of this parabola always lies on the x-axis, i.e., in the vertex we 
always get Qm — so that it will be interesting to have solutions attracted by it. Whether 
this parabola is open from above or below in the (x, Qm) plane depends on the sign of 
9(2o; a + <7 A ) 2 + 6(5c<ja + 2(x a ) which also determines the existence of xf so that if the critical 
points xf exist the parabola goes up and if they do not exist the parabola goes down. 
Since negative values of Qm are physically unreasonable, an open from below parabola 
does not represent a proper separatrix for the physically admissible region of the phase 
map. 

Close to the vertex of the separatrix, i.e., for: 

O [a: + 3(2^ A + q A + l)] 2 

Qm ~ [9(2^ + ^ + 6(5^ + 2^)]+^ (80) 

x = -3(2u x + a x + 1) + 5 X (81) 

with 5n M ,5 X <^1 the autonomous system becomes: 

dH 2uo x + a x + l 
~ 3 H 

dN 5 Qm 

dx a (S T 

~ 3(2w A + (7 A + l)- 



dVL M n 2cj a + a x + 1 



X2 



dN ~ 3(2lu x + a x ) 2 + 2(5cu A + 2a x ) 5, 



In the previous expressions, 8n M parametrizes the separation to the separatrix and is 
positive (negative) for points above (below) it, whereas 5 X gives the separation on the 
right (S x > 0) or on the left (5 X < 0) from the vertex of the parabola. From the equation 
for dH/dx we see that the trajectories will not be able to cross the vertex because dH/dN 
becomes singular at that point. On the other hand, it is easy to see from the equations for 
dx/dN and dVLu/dN that the vertex will always act as an attractor for the trajectories 
approaching from one side of the parabola so that, in the cases in which the parabola is 
open from above, there always exist trajectories that are attracted by the vertex of the 
separatrix. Whether the trajectories approaching the vertex are those going from above 
or below the parabola is determined by the sign of 2u x + a x + 1 as follows: 

i) 2u x + a\ + 1 > 0. In this case the trajectories approaching from below are at- 
tracted towards the vertex whereas those solution contained in the region above the 
separatrix are repelled by the vertex. 

ii) 2uj x + o x + 1 < 0. In this case, the solutions in the region above the separatrix are 
attracted towards the vertex whereas the trajectories below the parabola go away 
from the vertex. 
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To study the cases when the parabola is open from below we shall analyze the autonomous 
systems for x + 3(2c<j a + cr A + 1) = 5 X <C 1 and Qm <C 1- In that case we obtain 

(IX , , 5r 

3(2^ A + a A + l)- 



dN v Qm 

d£l]\j 

~~dW 



6(2u x + a x + l) (82) 



Thus, only if (2uj\ + cr A + 1) < the vertex of the separatrix is an attractor when the 
separatrix is contained in the region with Qm < 0. 

Finally, it remains to study the behavior of the equation of state as the trajectory ap- 
proaches the vertex of the separatrix. If we use the parametrization given in ( JHTl) again, 
we obtain that the equation of state becomes: 

2^ A + cr A + l 

w ~ 2 (83) 

<W 

Therefore, if the trajectory approaches the vertex from above the separatrix, the equation 
of state will evolve towards +oo (— oo) as long as 2o; A +cr A + l is positive (negative). On the 
contrary, when the trajectory goes to the vertex from below the separatrix the equation 
of state of the vector field goes to +oo (— oo) as long as 2u; A + a x + 1 is negative (positive). 

Then, if the parabola is open from above (3(2c<j a + (t a ) 2 + 2(5c<j a + 2<t a ) > 0), irrespectively 
of the sign of 2uj x + a x + 1 we have that the vertex acts as an attractor for some trajectories 
with solutions whose equation of state goes to — oo. On the other hand, when the parabola 
is open from below (3(2u; A + cr A ) 2 + 2(5c<j a + 2cr A ) < 0), only when 2u; A + a x + 1 is negative 
the vertex acts as an attractor and, in that case, the equation of state in the vertex goes 
to — oo. 



Finally, we would like to comment on the existence of certain models in which we can have 
critical points with Qm 7^ 0, 1. Those critical points can been found from (1751) by solving 
= with respect to x and for arbitrary values of Qm and, then, obtain the corresponding 
critical value for Qm from the equation = 0. In these critical points, one generally gets 
H = 0. The explicit expressions for these critical points are: 

x ± 2oo x - a x - 3 ± y/2(2cu A - a x - 3) [2(5a; A + 2a A ) + 3{2uj x + 

XQm ~ 4w A + 2a A + l 1 j 

o± 4 M H M + 3(2c A + a x + 1)] 

r A x n M ) 

Therefore, only models in which {2oj x — a x — 3) [2(5c<j a + 2a A ) + 3(2c<j a + cr A ) 2 ] is positive can 
contain these type of critical points. Moreover, this condition does not guarantee the existence 
of physically admissible values for Qm because one could, in principle, obtain both positive 
and negative values for Qm- However, these critical points cannot lead to accelerated solutions 
because the equation of state for the vector field in such points is identically zero, i.e., it behaves 
as a dust fluid. 
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To summarize the results of this section, we have shown that the features of the critical 
points for the case when we have matter in addition to the vector field remain unaffected, but 
the behavior of the separatrix presents novelties and, generally, in all the models we shall have 
attracting solutions with future singularities. Finally, we have shown that solutions in which 
Qm goes to some values different from and 1 are such that the equation of state of the vector 
field goes to zero, i.e., it asymptotically behaves as a matter fluid. 
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Figure 8: In this plot we show the regions where we have the vertex of the separatrix attracting 
trajectories with accelerated expansion (green shaded region) and where there are not solutions at- 
tracted by the vertex in which the expansion is accelerated (white region). We see that, in most of the 
parameter space, we have that the vertex attracts some trajectories in which the equation of state of 
the vector field diverges evolving towards — oo. The curve plotted in the graph separates those models 
in which the separatrix is open from above (outer region) and those in which the separatrix is open 
from below (inner region). Notice that these regions coincide with those in which the critical points 
P± exist (separatrix open from above) and they do not (separatrix open from below), as explained in 
the main text. 
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6 Local gravity tests 



The viability of any alternative theory of gravity is subject to its agreement with Solar System 
experiments, which provide very tight constraints on the so-called Parametrized-Post Newto- 
nian (PPN) parameters. These parameters are a set of quantities that characterize most of 
gravity theories at small scales and are extremely useful to measure deviations from GR. For 
the action (j2j), the PPN parameters are given by [18J: 



with: 



6 
A 



7 
a i 

0(2 
«3 



1 + AujA 2 (1 + ^±£) 
1 - AujA 2 (1 - 



1 + 



7(7 ~ 2) 
G 

4(1 - 7) [1 + 2eA] + 16uA 2 Aa 



^(3 + 7) + ^ 



3(1-7) 

e = Ci = 



4 ' 

1 + - £ A 



C2 — C3 — C4 



8uoA 2 Aa - 2- 



bA 2 



G 



56) 



(1 -AuuA 2 )(2e + a-2uj) 
(1 - AuA 2 )2e + 32u; 2 A 2 
1 

2A 2 a 2 - 2e[l - AA 2 {u + a)} 
2e(l - 37) + 2<t(1 - 2 7 ) 

(2cu + ( x) [(2 7 -l)(7 + l) + e(7-2)] 
-(2 7 -l) 2 (2cu + ( r + A) (1-2^+Aj A _^ 







A = 



Moreover, the effective Newton's constant is defined as: 



G 



eff 



G 



~1 -1 



-(7 + 1) + 6uM 2 ( 7 - 1) - 2AV(1 + 6) 



17) 



In the above expressions we have assumed G e // = 1 and A is the value of the vector field 
at Solar System scales (in units of AttG). The parameters (7,/?) are usually called the static 
PPN parameters and measure the space-curvature produced by a unit mass and the degree of 
nonlinearity relative to GR respectively. The parameter £ measures effects of preferred location 
whereas cti have to do with preferred frame effects. Finally, a.3 and Q are non-vanishing for 
theories in which the conservation of total momentum is violated. In GR, the PPN parameters 
are such that 7 — 1 = (3 — 1 = a\ = a,\ = a 3 = £ = £i = £2 = £ 3 = £ 4 = 0. On the 
other hand, for a general vector-tensor theory we see that there are neither preferred location 
effects nor violation of the total momentum conservations. However, these theories typically 
lead to preferred frame effects (as expected because of the presence of a vector field) as well as 
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deviations from GR for the static PPN parameters. Current observational limits on the PPN 
parameters impose very stringent limits on modified gravity theories because they do not allow 
much deviation from GR, i.e., GR agrees with local gravity tests with very good precision [18J: 

7-1 < 2.3 x 1(T 5 

0-1 < 2.3 x 1(T 4 

ax < lO" 4 

a 2 < lO" 4 (10- 7 ) (88) 

In order to obtain constraints on the vector field from these limits we linearize the PPN pa- 
rameters given in (|86|) as follows: 

7 - 1 ~ — [2(e -u)+cr] A 2 
e 

. 1 (2u;-a)(2e-2u; + a)(4e-2u; + <7 A ) 2 

0—1 ~ — ; n -A 

4e 2 

16cu(2e + a) „ 2 
ai ~ -A 

f [(2L + ( T) 2 (I-f)-4o-]A 2 A^O 
a 2 ~ ^ (89) 
[ *f [2(e + w) + a] A 2 A = 

Therefore, we can set that, typically, the vector field at the Solar System scale will be con- 
strained to be A < 10~ 2 , for models in which all the parameters are order unity. Let's remark 
that this value of the vector field does not need to coincide with its cosmological value. 

On the other hand, the linearized Newton's constant is given by: 

\2uo-a) 2 



G ~ 1 



+ 4(<7 - u) 



A 2 (90) 



Then, if we use the existent limits on its time-variation G/G < 10~ 13 yr -1 together with the 
constraints on the vector field obtained above, we can also set bounds on the cosmological time 
variation of the vector. 

Although a general vector-tensor theory will be constrained by the aforementioned limits, 
there is a number of models whose parameters satisfy certain relations so that some of the PPN 
parameters could be identical to those of GR and, thus, pass the limits given in (|88|) . Indeed, 
in [7] is shown that there exist a total of 6 models which are indistinguishable from GR by 
means of local gravity tests for any value of the background vector field, namely: {a = 0, e, A}, 
{a = — 4A = — 4e}, {a = — 3A = — 2e} and {cr = me, with m = 0, —2, —4}. Moreover, in that 
work, a detailed analysis of the stability of these models is performed. 
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7 Conclusions and discussion 



In this work we have developed a general study of the cosmological evolution of a vector field 
non-minimally coupled to gravity and without potential terms. We have given the evolution of 
this vector field in terms of the parameters of the theory for the different epochs of the expansion 
history of the universe, namely: inflation, radiation dominated era and matter dominated era. 
We have shown that it is possible to obtain a wide variety of behaviors for the evolution of the 
vector field by suitable choices of the parameters. In particular, we have obtained conditions 
for the parameters so that the vector field energy density grows or decays with respect to that 
of the dominant component. Moreover, conditions to have scaling evolution have also been 
calculated. 

The case of a universe dominated by the temporal component of the vector field has been 
studied in detail. We have obtained an autonomous system describing the evolution of the 
Hubble expansion rate and the vector field. The general features of the phase map have been 
given and all the critical points have been appropriately characterized. For those points that 
act as attractors we have obtained the general conditions under which the vector field gives rise 
to accelerated expansion. 

To study the viability of these theories as dark energy candidates we have performed an 
analysis of the field equations together with Einstein's equations when the universe contains 
matter in addition to the vector field. Then, we have identified solutions which allow a transition 
from matter-domination to vector-domination with accelerated expansion. We have also shown 
that these models generally contain future singularities in which the equation of state diverges 
and that most of the models can give rise to periods of accelerated expansion. 

In addition to the general results commented above, we have also found a wide variety of 
particular model examples with interesting properties, thus: models with late-time attractors 
with Qm ^ 0, 1 all with equations of state for the vector field resembling that of non-relativistic 
matter. This type of models could provide vector dark matter candidates. We have also found 
models with early time accelerated solutions which are unstable thus giving rise to possible finite 
inflationary periods. On the other hand, the dark energy model proposed in [12] corresponds 
to the parameters e = —1/4, A = —1/2, o = 1, u = 0. Concerning the temporal component, 
we have that a\ = —2 and oj\ = and, if we look at Fig. [TJ we see that it lies in the 
line of scaling behavior during radiation (what allows to solve the coincidence problem in a 
natural manner) and in the growing region for the matter era, in agreement with the result 
found in [12] . Moreover, it belongs to the region I of the classification summarized in Table 
CD for the vector dominance case so that Pq and P + behave as attractors whereas P_ and the 
separatrix repel the trajectories. The singularity found in [12] corresponds to the trajectories 
approaching the separatrix in the case when matter is present. This is also in agreement with 
the results of this paper because we have that the separatrix for this particular model is given 
by VLm = j2(x — 3) 2 so that it is open from above. Moreover, we have that 2^ + 0^ + 1 = — 1 < 
so that, according to the discussion about the separatrix in section 5, the trajectories above 
the separatrix will be attracted by its vertex. Then, as the vector field evolves as oc t a with 
a = i(-3 + \/3) in the matter dominated era, we have that the initial conditions are given 
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by VLm — 1 and x ~ § a ~ 0.69 so that the cosmological solution corresponds to one of 
the trajectories approaching the vertex of the separatrix and, therefore, leading to a future 
singularity. Finally, the electromagnetic dark energy model proposed in [13] has parameters 
uj = a = 0, A = 1/6, e = — 1/4 and coincides with the only model whose temporal component 
has constant energy density for all the cosmological epochs (see Fig. [T]). Moreover, if we look 
at Fig. [3] we see that it belongs to Case Ilia in which both critical points P + and P_ coalesce 
and they coincide with the sepratrix. In [13] it is shown that this model, indeed, behaves as a 
cosmological constant irrespectively of the Hubble expansion rate so that it leads to a de-Sitter 
universe once it dominates over the matter component. 

In the last section of the paper we have given the expressions for the PPN parameters of the 
general vector-tensor theory and used current limits on the PPN parameters to constrain the 
value of the vector field at small (Solar System) scales, which has been found to be A Q < 1CT 2 . 
However, one cannot, in principle, use this bound to put limits on the cosmological value of 
the vector field because both values do not need to coincide, unless the power spectrum of the 
vector field happens to be scale- invariant. In the general case, one would need to know the 
mechanism that originated the primordial power spectrum of the vector field so that we can 
know the expectation value of the vector field at Solar System scales and, thus, compare to the 
experimental limits. 

We would like to remark that, throughout this work, we have focused on the viability 
of vector-tensor theories as candidates for dark energy just for the first stage, i.e., we have 
given the general conditions under which the homogeneous part of a vector-tensor theory can 
lead to late-time acceleration. However, to propose a serious candidate one should also study 
the perturbations of the corresponding model and check the presence of instabilities both at 
classical and quantum levels. The absence of classical instabilities is required in order not 
to have exponentially growing modes that could spoil the predictions of the model for the 
zero-mode. The issue of quantum instabilities has to do with the existence of modes with 
negative energy (ghosts) so that non-linear interactions of the field might produce an unlimited 
number of such particles. The existence of such instabilities has already been studied in [7] for 
those vector-tensor theories which are indistinguishable from GR at small scales by considering 
perturbations in both the vector field and the metric. In [23], it was shown the existence of 
instabilities for the longitudinal component of the vector field in some particular cases. In [5], 
a detailed treatment on the stability of the vector field perturbations for general vector-tensor 
theories with a potential term was performed. However, a complete study by considering metric 
perturbations in addition to vector perturbations remains to be done, although such a study 
is far from straightforward because of the large number de cases involved. In any case, as 
suggested in pi] where general N-forms are studied, the existence of singular classical solutions 
is related to the presence of ghosts so that one would expect that most of the vector-tensor 
models present instabilities (at least at the quantum level) because, as shown in this work, the 
existence of singular solutions at the classical level is a common feature in these models. 

As a final comment, although throughout this work we have focused on vector fields as 
candidates for dark energy, we would like to mention that, because of their generality, the 
results given in the paper can also be used in other cosmological contexts in which vector fields 
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could play a relevant role. 
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